Polaron formation: Ehrenfest dynamics vs. exact results 
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We use a 1-dimensional tight binding model with an impurity site character- 
ized by electron-vibration coupling, to describe electron transfer and localiza- 
tion at zero temperature, aiming to examine the process of polaron formation in 
this system. In particular we focus on comparing a semiclassical approach that 
describes nuclear motion in this many vibronic-states system on the Ehrenfest 
dynamics level to a numerically exact fully quantum calculation based on the 
Bonca-Trugman method [J. Bonca and S. A. Trugman, Phys. Rev. Lett. 75, 
2566 (1995)]. In both approaches, thermal relaxation in the nuclear subspace is 
implemented in equivalent approximate ways: In the Ehrenfest calculation the 
uncoupled (to the electronic subsystem) motion of the classical (harmonic) os- 
cillator is simply damped as would be implied by coupling to a markovian zero 
temperature bath. In the quantum calculation, thermal relaxation is imple- 
mented by augmenting the Liouville equation for the oscillator density matrix 
with kinetic terms that account for the same relaxation. In both cases we 
calculate the probability to trap the electron in a polaron cage and the prob- 
ability that it escapes to infinity. Comparing these calculations, we find that 
while both result in similar long time yields for these processes, the Ehrenfest- 
dynamics based calculation fails to account for the correct timescale for the 
polaron formation. This failure results, as usual, from the fact that at the 
early stage of polaron formation the classical nuclear dynamics takes place on 
an unphysical average potential surface that reflects the otherwise-distributed 
electronic population in the system, while the quantum calculation accounts 
fully for correlations between the electronic and vibrational subsystems. 
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I. INTRODUCTION 

Electron transfer (ET) between molecular systems has long been recognized as a key process 
in many research fields of chemistry, physics and biology 1-5 . Many of its aspects are described 
by the Marcus theory 6 , which has been extended to describe such areas as artificial solar-energy 
conversion 7-9 and molecular electronics 10,11 . 

The Marcus theory relies in an essential way on electron-vibration interaction. The initial and 
final states of the electron transfer process are fully equilibrated polarons localized on different 
sites, and transitions between them is evaluated within the assumptions of transition state theory. 
Motion in an extended system is assumed to be a succession of hopping steps, each described as a 
Marcus process. In the other extreme limit, electronic motion in a frozen lattice, the electron moves 
within its energy band, most simply described using a tight binding model. In between these limits, 
electron-phonon interaction and band motion can change the electron's character from being weakly 
perturbed by electron-phonon scattering to polaronic motion whose discrete representation is the 
succession of hopping processes described above. 

In the present paper we are interested in situations where electronic band motion competes on the 
same timescale with polaron formation, so that the dynamics of the latter process has to be considered 
explicitly. Such considerations are relevant to recently studied models of photovoltaic cells 12,13 , where 
electrons (or holes) are injected at some location in the system and a useful process is defined by their 
absorption at another (e.g. an electrode surface). The yield of such processes, determined by the 
competition between electronic motion and loss processes 13 (e.g. carrier recombination) is expected 
to be sensitive to electron-phonon interactions, and in particular to transient polaron formation. 

Exact treatment of such coupled many-body systems is difficult, and it is tempting to resort to 
approximations such as the semiclassical mean field (Ehrenfest) dynamics. In this approximation 
the electronic wavefunction \&(r, t) (r represents the electronic coordinates and t is the time) evolves 
under a time-dependent Hamiltonian defined by a classical nuclear trajectory, schematically repre- 
sented by a nuclear coordinate R(t), while the latter is obtained by solving the Newton equation 
for the nuclear motion with a potential in which the vibronic coupling V(r, R) is replaced by its 
instantaneous expectation value V(R, t) = {9(r,t)\V(r,R)\^(r,t)). Such an approximation, essen- 
tially a dynamical extension of the Born-Oppenheimer (BO) approximation, is expected to perform 
well when the electronic motion is fast, throughout the relevant electronic subspace, relative to the 
nuclear dynamics. Its failure in describing processes in which transitions between BO electronic 
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adiabatic states take place on the same timescale as nuclear motions is also well known. Indeed, 
in the analogous case of electron solvation in polar liquids such non-adiabatic processes have been 
addressed with the necessary accounting for the quantum nature of the nuclear motion 14-16 , usually 
within the surface hopping methodology 15-17 . Still, because Ehrenfest dynamics is so easy to im- 
plement and to use, it is of interest to assess its performance as an approximation to exact results 
in the context described above 18 . This is the purpose of the present paper. Using a model that is 
simple enough to solve up to any desired level of accuracy, we focus on two observables: the extent 
of the polaronic localization and its formation time, and compare results obtained from the semi- 
classical Ehrenfest dynamics approach to the exact, fully quantum, results. For model parameters 
that support polaron formation we find that, while the Ehrenfest calculation yields a similar final 
state as the exact one, it predicts a polaron formation time that is an order of magnitude longer 
than the exact result. This implies that Ehrenfest dynamics cannot be used as a reliable tool for 
assessing polaronic effects in such systems. This does not exclude its possible applicability in larger 
systems with higher temperatures with many more nuclear degrees of freedom, but indicates that 
its use should be exercised with caution and after performing suitable benchmark calculations. 

We start by formulating the basic Hamiltonian model and comparing the different predictions of 
the quantum and the semiclassical descriptions in Section II. In section III we define the population 
operator and the population formation time. Numerical calculation and discussions are given in 
sections IV and V, and a conclusion follows. 

II. THEORETICAL MODEL 

We consider an n+l-site tight-binding free electron model (below we take n=4) coupled to a 
system of Harmonic oscillators. We assume that only one oscillator (henceforth referred to as the 
"primary" vibration) directly couple to the electronic systems. The others ( "secondary" phonons) 
constitute a thermal bath that affects relaxation in the primary system. The Hamiltonian of the 



FIG. 1: The model chain including 5 sites. The electron-vibration interaction occurs on the middle site 2. 



whole system is 



H = H s + H b + Hsb , (1) 

4 3 

H s = £/cJq + V }X4 C W + 4+i c i) + ^o4,d + a 2 4c 2 (4 + rf o) , (2) 
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H B = ^ hjj s d\d s , (3) 

s=l 
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^ = ^A s (44 + 44) • (4) 
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Here ifs corresponds to the electronic system (described by the creation and annihilation oper- 
ators for each site /, cj, together with the primary vibration, of frequency uq described by the 
creation and annihilation operators rfj, <^o- H B describes the secondary phonon bath (d\, d s are the 
creation and annihilation operators for phonon of frequency u s ) and H$b is the coupling between the 
primary vibration and secondary phonons. E\ is the on-site energy level of site /, V is the coupling 
parameter associated with electron tunneling between nearest neighbor sites. The parameters a 2 
and X s correspond to the coupling between the electronic state at site 2 and the primary vibration, 
and between the primary vibration and secondary phonons, respectively. 



A. The quantum approach 

In the quantum approach, the evolution described by the Hamiltonian (2) is treated essentially 
exactly 18 , using the basis set {|n, v >} where n and v denote the electronic state localized on site n 
and the vibrational state v of the primary oscillator. In the numerical calculation we truncate the 
set {u} at some value, v max and test for convergence as v max increases. 



The coupling to the thermal bath is treated in the master equation approach: The density matrix 
pT of the whole system is assumed to keep the form px = Ps®Pb, where ps, the density matrix of the 
thermal bath (secondary phonons), is assumed to remain in equilibrium at the ambient temperature. 
Then the quantum master equation (QME) for ps is 

iU ^W^ = i H SiPs(t)}-ifrlol4d ps(t)+p s (t)4d -2d oPs (t)4}/2 , (5) 

where 70 is the vibrational relaxation rate induced by the vibration-phonon bath coupling. It is 
equal to the imaginary part of the vibration self energy E given by 

7o(w)/2 = ^Im{S vibration (w)} = - ^ \K\ 2 S(hu - hu s ) . (6) 

s 

In the basis chosen, Eq. 5 takes the form 

ih dPnV '^' f = [H S , p S (t)]nv,n'v> - ^7o0 + v') Pnv,n'v'^n,n' /2 + ikjoy/v + Wv' + lp n ,v+l,n' ,v> +1 , (7) 

where the first term in the right side of Eq. 7 comes from the contribution of primary system 
Hamiltonian Hs, and second part comes from the first two terms of the bracket on the right side of 
Eq. 5, and the last term involves energy transfer from the higher to the lower vibrational levels 31 . 
This equation will be solved numerically. 



B. The semi-Classical approximation 

The dimensionless displacement of the single primary oscillator is approximated in the semiclassical 
approximation by a time- dependent configuration q(t) =< 4(t) + d (t) > as 19 ' 20 

1 /** 

q(t) = - J drD r {t-r)a 2 <c\{r)c 2 {r)> , (8) 

arising from the interaction among the electron, the active vibration and the phonon bath; here D r 
is the retarded green function of the active vibration. 

Equation 8 assumes that the displacement coordinate q responds to the average electron popula- 
tion (in the present calculation, on site 2). This is a mean field description akin to the Ehrenfest 
approximation, that is to be tested in the calculations described below. 



Using the wide-band approximation 

D r (oo) = ; (9) 

w + a;o + ^7o/2 u - cu + ijo/2 

and its Fourier transform 

D r (t) = l [ e ^-^l 2 ) t - e (-^o-7o/2)t] = ^2stn(u t)e-^ ot/2 , (10) 

and substituting Eq. 10 into Eq. 8, we get 

qit) = -\ [ drsin[u Q {t - r)]e- 7o(t - r)/2 a 2 < 4(r)c 2 (r) > . (11) 
"■Jo 

Here 70 has been defined in Eq. 6, and neglecting the real part of the vibration self energy gives 
the solution for q(t). 

Finally replacing + d in the electron- vibration coupling Eq. 2 by q(t) (Eq. 11), we get the 
effective semiclassical electronic system Hamiltonian as 19 ' 21 ' 22 

4 3 

H eff = Yl £lC i Cl + V J2(4ci+i + 4+iCi) + F{t)c\c 2 , (12) 

Z=0 1=0 

with 32 

F(t) = a 2 q(t) = -H^l f dTsin[uj (t - r)]e^ o( '" r)/2 < 4(r)c 2 (r) > . (13) 
" Jo 

The system density matrix p$, which in this semiclassical approximation is derived from the 
Hamiltonian H e ff in Eq. 12, can be solved using the Liouville equation 

ih^ = [H eff , Ps }. (14) 

III. POPULATION DISTRIBUTION, POPULATION FORMATION TIME AND 
ELECTRON- VIBRATION COUPLING ENERGY 

The on-site electronic population at any time t is 

4 

p,(t) =< 4(t) Cl (t) > ,with J2 P ^ = 1 • ( 15 ) 

1=0 

Since the time-dependent values P; oscillate, it is better to show these values using a coarse grained 



time-dependent average value 

-I pt+AT 
H—AT 



Pi(t) = 2^ / <ItPi{t) . (16) 



Below we use AT = 50 fs. In the following we will use p as the time-dependent average values for 
the populations. 

The population formation time of P2 can be defined as the time point at which population P 2 
reaches a certain value. Experimentally one can define the "formation time" as the time at which 
the target population reaches ~ (1 — e _1 ~ 0.76) of its final value 23 . Thus we use the criterion 

P 2 (r p ) = P 2 °°(l-e- 1 ). (17) 

P 2 °° is the time-averaged value of P 2 m the long time limit and t p is the population formation time. 
The electron-vibration coupling Ep are 



E p 



«9 < c\c2{dl + do) >, quantum method, 

(18) 

F(t) < cJ,c 2 >, semiclassical method. 

IV. NUMERICAL CALCULATION 

For the initial numerical simulation, we set E\ = for / = 0, 1, 3, 4 with e 2 = —0.2eV(e2 is lower than 
the other site energies), V = O.leV, v max = 9, a 2 = 0.0707eV, u = O.leV, 70 = 0.04e^. We will 
vary those parameters to examine more cases. All the results described below use the initial condition 
c li c n = ^n,o f° r the electronic state, and v — (ground vibrational level) for the primary phonon. For 
our purpose-comparing the quantum calculation and the mean-field semiclassical approximation- it 
is sufficient to consider the zero temperature case. 

For the quantum method, we need to set v max large enough to assure convergence of the calculation. 
Because we are at zero temperature, the energy transfer only happens from the higher levels to its 
nearest lower level and a smaller v max = 3 is enough (in the supporting information Fig. SI shows 
that when v max is larger than 2 the average population P 2 is similar). For non-zero temperature, 
there will be transfer from lower to higher levels, and larger v max will be needed. 
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V. POPULATION LOCALIZATION AND POLARON FORMATION 

The electron- vibration coupling is allowed only at site 2 and the system is coupled to the bath by 
a single primary vibration. As shown in Fig. 2, for a small V, much of the electron population will 
localize on site 2, forming a local polaron with larger population. The quantum results (Section II) 
are shown on left panel and the semiclassical results (section III) are displayed in the right panel. The 
respective results qualitatively similar:^ (polaron population) rises and saturates at values that are 
similar in both approaches. However the rise time obtained from the quantum calculation is much 
shorter than its semiclassical counterpart: P 2 reaches its maximum value at 2 ps (picosecond) in the 
left panel, while in the right panel it takes about 15 ps. This difference demonstrates the shortcoming 
of the semiclassical approximation, or rather-its reliance on the mean field approximation: The 
localizing phonon moves on a potential surface that is substantially different from what it actually 
experiences once the localization process has started. 

A. Influence of the tunneling amplitude V 

The coupling V between nearest neighbor sites is important for population localization. For small 

V the coherent transfer between the different sites is weak and the population tends to localize, 
forming a polaron on site 2 which is coupled to the primary vibration and has been given a lower 
energy. For large V the population tends to equalize on neighboring sites, showing an oscillatory 
behavior (here only the average is shown ). As shown in Fig. 3, the average population P 2 on the 
impurity site decreases with increasing V, although the polaron is in principle formed on site 2. For 

V = l.OeV, the banding or delocalisation energy dominates over the electron phonon and impurity 
energies. In this limit, the charge population will tend to reach a uniform average value. Recall that 
the original small-polaron model of Holstein 24-27 was developed for narrow-band (small V) materials. 

B. Effect of the electron-vibration coupling 

a 2 is the coupling strength between the electronic motion and the vibration at site 2. In the 
quantum method, electrons can evolve from the other sites to the vibronic states through this 
coupling as shown in the last term on the right side of Eq. 2. With larger a 2 , more population will 
transfer and form a polaron with large population on the middle site. For the semiclassical case, the 
same thing will occur through Eq. 13, and some population will be localized to form a polaron. As 
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shown in Fig. 4, the population P 2 in the steady state increases with this coupling strength. The 
more rapid population relaxation obtained in the quantum treatment occurs because in the quantum 
analysis, the actual population at site 2 is used, while in the semiclassical treatment, it is necessarily 
only the average population which drives the coupling. 

In Figs. 5(a) and 5(b), P 2 °° (defined in Eq. 17) is plotted in 3-D by changing the coupling V 
and the electron-phonon coupling a 2 . With a small V, P 2 °° is almost around 0.5 33 . Choosing larger 
values of V, decreases the site 2 population as can be expected, as shown in the corners of Figs. 5(a) 
and 5(b). However increasing a^, the population builds up again, and we can say that the polaron 
population value is mainly determined by the parameters V and a^- 

C. Comparison of population formation times 

In Figs. 6(a) and 6(b) the population formation time is shown in 3-D by changing both V and 
a 2 • The formation time obtained by the semiclassical method is longer than the quantum method. 
The population formation time decreases with electron-phonon coupling a 2 but then saturates. It 
also decreases with coupling V, but there is a turnover in the formation time beyond a V~0.1eV. 
This is due to the fact that for large V, the excitation relaxes into a delocalized population which is 
no longer strictly speaking a localized polaron. The population distribution is roughly constant in 
this limit as shown in Fig. 3. The formation time should now be referred to simply as population 
relaxation time. 

D. The short time dynamic of polaron formation 

In Fig. 7(a) we compare the time- dependent dynamic processes for P 2 and Ep (Eq. 18). Both of 
them reach their steady state very quickly. When P 2 reaches its maximum, E p is at minimum, and 
vice versa. This represents the damping of a coherent oscillator where displacement and population 
keep their phase delay throughout. The picture (Fig. 7(b)) is somewhat different in the semiclassical 
case, P2 and its F2 are also delayed but we now see beats due to the interference of waves scattering 
from an oscillating potential. Both quantities reach their steady state more slowly in this case. 

VI. CONCLUSION 

In this paper we have compared a fully quantum and semiclassical model. The latter based on 
Ehrenfest dynamics for the calculation of the polaron formation process. We used a 1-dimensional 
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tight binding model that includes electron-phonon coupling but only one site, the so called polaron 
trap, the "trap phonon"is in turn coupled to a thermal bath which allows the system to relax into 
"an equilibrium" . 

The results from both methods show qualitatively similar behavior, however with markedly dif- 
ferent timescales: the population formation time obtained from the quantum calculation can be 
10 times faster compared to the semiclassical method. This discrepancy becomes smaller with in- 
creasing intersite coupling V (for large V no localized polaron is formed) since the classical limit is 
reached for V ^> oj . 

The different relaxation times obtained in the quantum and the semiclassical calculations result 
from the use of mean field approximation in the latter. In this approximation, the primary oscillator 
responds to the average occupation of site 2 which effectively makes it move on a potential surface 
that is markedly different (less binding) than the one it experiences once localization is initiated. 
Localization on this average potential is slower. This averaging assumption may be justified when 
the electron motion is much faster than the vibration, that is, for large V (V ^> faujo). Another minor 
difference between the two calculations is the use of a master equation in the quantum calculation, 
and an essentially equivalent Langevin equation in the semiclassical one, to describe the relaxation 
of the primary phonon. These relaxation schemes are equivalent, provided that care is taken to use 
parameters that imply the same relaxation rate in both cases, as was done here. 

The conclusion would also apply to the "exact" simulations of Kopidakis et al 28 and the many 
other papers where the vibrations are also treated with semiclassical dynamics. It would seem that 
depending on electron bandwidth, the formation time is considerably underestimated in these works, 
perhaps by an order of magnitude or more. The semiclassical results are essentially in agreement with 
the conclusions reached by Emin and Kriman 29 ' 30 using the Holstein diatomic polaron lattice model. 
These authors showed that population localization and polaron formation depend critically on the 
ratio of the tunneling energy V and the width of the Bloch phonon dispersion, which effectively 
plays the role of the dissipation term since outward traveling phonon Bloch waves will not return. 
One final but very important point. We have shown that for simple relaxation models, the true 
thermodynamic equilibrium is not necessarily the one reached in the steady state. The state reached 
with arbitrary start conditions can, in terms of energy, be a metastable state the populartion on site 
2 (lower energy ) never exceeds however small we make V. Introducing even a small coupling on the 
other sites then allows the system to reach the true ground state and now the population can climb 
up to 1. This very interesting point needs to be investigated in detail and in particular at finite 
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temperature. 
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FIG. 2: Average population distribution on different sites I (I = 0, 1, 2, 3, 4) shown as a function of time, ei = (I = 0, 1, 3, 4), 
e 2 = -0.2eF, V = O.leV, aj = (I = 0, 1, 3, 4), a 2 = 0.0707eV, &j = 0.W, fryo = 0.04eV. Due to symmetry P =Pa; Pi=Ps- 
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FIG. 3: Average population P2 as a function of time with different nearest neighbor tunneling parameter V. ei = (I = 0, 1, 3, 4), 
£2 = -0.2eV, Q; = (I = 0, 1, 3, 4), a 2 = 0.07076^, froj = O.leV, ft 7o = 0.04eV. 
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FIG. 4: Average population P2 shown as a function of time with different electron- vibration coupling parameter <»2- £1 = 
{I = 0, 1, 3, 4), e 2 = -0.2eV, V = O.leV, a t = (I = 0, 1, 3, 4), &j = O.leV, h-y = 0.04eV. 
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Quantum method 




(a)quantum method 

Semiclassical method 




(b)scmiclassical method 

FIG. 5: P 2 °° (population on site 2 in the steady state) shown as a function of the nearest neighbor site coupling parameter V 
and electron-phonon coupling o>2. £i = (I — 0, 1,3,4), £2 — — 0.2eV, hojo = O.leV, /ryo = 0.04eV. Quantum method used for 
panel (a) and semiclassical method used for (b). 



17 




(a)quantum method 




(b)scmiclassical method 

FIG. 6: Population formatiom time for P2 shown as a function of the nearest neighbor site coupling parameter V and electron- 
phonon coupling oti. ei = (I — 0, 1,3,4), £2 = —0.2eV, fkoo = O.leV, hyo = 0.04eV. Quantum method used for panel (a) and 
semiclassical method used for (b). 



Quantum method 
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Time [fs] 

(b)scmiclassical method 

FIG. 7: (a): Population P2 and electron- vibration coupling energy Ep are showing as a function of time; (b): Population P2 
and the electron-phonon coupling energy F2(t) are showing as a function of time, ei = (I — 0, 1,3,4), £2 = —0.2eV, Q; = 
(I = 0, 1, 3, 4), «2 = 0.0707ey hwo = O.leV, = 0.04eV '. Quantum method used for panel (a) and semiclassical method used 
for (b). 



